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Abstract 

Graeffe iteration was the choice algorithm for solving univariate polynomials 
in the XlX-th and early XX-th century. In this paper, a new variation of Graeffe 
iteration is given, suitable to IEEE floating-point arithmetics of modern digital 
computers. 

We prove that under a certain generic assumption the proposed algorithm con- 
verges. We also estimate the error after A'^ iterations and the running cost. 

The main ideas from which this algorithm is built are: classical Graeffe iter- 
ation and Newton Diagrams, changes of scale (renormalization), and replacement 
of a difference technique by a differentiation one. 

The algorithm was implemented successfully and a number of numerical ex- 
periments are displayed. 
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1 Introduction 

Many present day numerical algorithms have originated in highly ac- 
claimed methods dating from last century or even earlier. Such is certainly the case of 
Euler's method or of Newton's method, whose numerical and theoretical consequences 
still impact us today [15, 33, 34, 35, 36, 37, 39]. 

Graeffe's classical method for finding simultaneously all roots of a polynomial was 
introduced independently by Graeffe, Dandelin and Loba- 
tchevsky [11]. Its simplicity, as well as importance throughout last century indicate 
its potential as an effective numerical algorithm. 

Surprisingly, Graeffe's method has not received much attention in 
present day numerical computations. Very few modem discussions about it or its ap- 
plications can be found. See the review by V. Pan [28], and also [2, 5, 6, 8, 16, 21, 22, 
24, 27, 29, 32]. 

One of the main reasons for Graeffe's lack of popularity stems from the fact that 
its traditional form leads to exponents that easily exceed the maximum allowed by 
floating-point arithmetic. Other reasons, such as the "chaotic" behavior of the argu- 
ments of the roots of the iterates contribute to such stigma. 

Also, Graeffe iteration is a many-to-one map. It can map well-conditioned polyno- 
mials into ill-conditioned ones, as pointed out by Wilkinson in [40]. We shall refer to 
this as 'Wilkinson's Deterioration of Condition'. 

In this work we present a version of Graeffe's algorithm, which is well suited for 
floating-point arithmetic computations. Furthermore, it has exceUent complexity and 
memory allocation characteristics. Our method computes both the moduli and the 
argument of all the roots, provided that certain generic conditions are satisfied. These 
claims are backed by our theoretical results presented in the next section, and proved 
throughout the paper, as well as the numerical experiments presented in the end of the 
paper. 

The main ingredients in our approach are the foUowing: 
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• The idea of renormalizing the relevant operations at each iteration step, akin to 
what is done in dynamical systems and physics [19, 23]. 

• The idea of using the differential of our RenormaUzed Graeffe iteration as a way 
of keeping the information concerning the argument of the roots. This will allow 
us to avoid the harmful effects of Wilkinson's 'deterioration of condition', as 
discussed in Section 4.6. 

• A renormalized version of Newton's diagram that allows us to recognize and 
locate pairs of conjugate roots, as well as roots of higher multiplicity. 

The first idea mentioned above was developed in our earher work [22], which in a 
certain sense laid the conceptual framework for our present approach. It is not however 
essential in understanding the proofs presented herein. 

The second ingredient mentioned above is explained and motivated in Section 4. 
In rather vague terms it could be compared to the advantage of using derivatives, when 
those are available, as compared to using differences. This idea can be traced back to 
Brodetsky and Smeal [4] in 1924, in a more ad-hoc fashion. We are not aware of recent 
applications of that method in modem Uterature. 

Finally, the concept of Newton's diagram, as well as the power of Graeffe's method 
was present throughout Ostrowski's masterpiece [25]. While writing the present paper 
we could not help but wonder what would have been the outcome of that research if he 
had available at that point the present day technology of high speed computers. 

We wish to thank two anonymous referees for their cormnents and for suggesting 
some extra references such as [4], which we were not aware of in the first draft. 



1.1 IVIain Result 

We will introduce an algorithm for solving real and complex univariate polynomials. 
The following genericity condition will be required at input: 

Definition 1. A real polynomial / will be called circle free if, and only if, for any 
couple (, ^ of distinct roots of /, one has either |CI 7^ 1'?!^ or ^ = ^. 

Definition 2. A complex polynomial / will be called circle free if, and only if, for 
any couple C. C of distinct roots of / one has \C,\ ^ |^| . 

It is obvious that given any real polynomial /, one can obtain a circle free polyno- 
mial by (pre)composing it with a conformal transform of the form: 

<yj: C ^ C 

X ^ ipix) = '^co^^-si"^ 

and then clearing denonainators; for all but a finite number of ^ € (— tt, tt], the resulting 
polynomial 



f{x) = {xsme + coseff 



X cos f — sm ( 
a; sin 4- cos ( 
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is circle free, where d is the degree of /. 

Tangent Graeffe Iteration will be shown to converge for all circle free polynomials; 
Given an arbitrary polynomial, one can first find all zero roots (in the obvious way), ap- 
ply a random conformal transform, then Tangent Graeffe Iterations, and finally recover 
the roots of the original polynomial. 

When counted with multipUcity, the roots of a circle free polynomial can be canon- 
icaUy ordered by: 



If we assume that all the arithmetical operations are performed exactly (including 
transcendental), the mathematical properties of the algorithm can be summarized by: 

Theorem 1. Let f be a real ( resp. complex) circle free degree d polynomial, not 
vanishing at 0. Denote by Q the vector of all the roots of f with multiplicity canonically 

ordered as above. 



Then, a total of N iterations of Renormalized Tangent Graeffe (Algorithm 6) pro- 
duces C^^^ e C'', such that 



The running time for each iteration is 0{(P ) exact arithmetic operations ( including 
transcendental operations). The relative truncation error bound in each coordinate 
after N iterations is , where C depends on f. 

1.2 What the Graeffe Iteration is; Its historical weaknesses 

In this section we shall briefly review the main ideas behind the method and describe 

also some of its weaknesses. 

Graeffe iteration maps a degree d polynomial f{x) into the degree d polynomial 



If Cii C2j • ■ • Cd are the roots of /, then the roots of Gf are Ci i Cl) • • • Cd 

Assume that g = G^f is the A'^-th iterate of /. Then, assuming that / is monic, 
the coefficients of g{x) = go+ gix H h gd.x'^ satisfy: 




2.1. If|Ci| = |Ci+i|thenCi = Ci+i- 

2.2. If i = 1 or |Ci-i| < ICil then Im Ci > 0. 



— ^ ^ (asN ^ oo). 




9o = {-^y-'^d 



9d = o-Q 



= 1 
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where aj is the j-th elementary symmetric function. In the particular case that |Ci| < 
IC2I < • • • < |Cd|> we can further approximate 

.90 = i-infcf-.-cf 

9. - (-i)'^-^cr...cr 

9d = (To = 1 
Hence, it is possible to determine 

^2^^ ^ 9j 
^ 9j+i ' 

We stress two main weaknesses. The first big weakness of classical Graeffe itera- 
tion is coefficient growth. As the coefficients of gj grow doubly exponentially in the 
number of iterations, the exponent (not the mantissa) of the floating-point system gets 
overflowed: 

Example 1. Let / have roots 1, 2, 3, 4. Then the A'^-th Graeffe iterate of / has roots 
1, 2^ ,3^ ,4^ . The coefficient go is 24^ . If iV = 8, then go is approximately 
1.68 X 2^^^^, while IEEE double precision numbers used in most modern computers 
cannot contain floating point values more than 2^"^^, since the exponent is represented 
by 11 bits (sign included) [10]. (As a matter of fact, the representation is a little more 
complicated, as it allows for 'subnormal' numbers [7]). Therefore, we would have an 
overflow when computing the 8-th Graeffe iterate of /. 

Example 2. On the example above, assume that / would have an additional root 1.01. 
Namely, 

f{x) = (.T-l)(a;-1.01)(a;-2)(a;-3)(.T-4) 

= -24.24 + 74.5.T - 85.35.t2 + 45. la;^ - ll.Qla;'' + 

We will show that 8 Graeffe iterations are not enough to compute the first root (namely 
1) to an accuracy of 10~*. However, as shown in Example 1, 8 iterates are enough to 
overflow the IEEE double precision number system. 
Indeed, the first root is obtained as: 

1^2^ SO _ 1.01'' 24^ 

^ - 91 (12™ +1.012" )242~ 

~ 1 + 1.01-2" 

~ 0.927 (for AT = 8). 



Thus, 



C~ 1-2.9 X 10"^ 



The error obtained is therefore larger than 10 ^. 
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The introduction of the idea of renormalization allows us to avoid coefficient 
growth, and replace a diverging algorithm by a convergent one (See Section 2 and 
also [22]). Alternative approaches for the number range growth are suggested in [9] 
and in [8]. 

A certain geometrical invariant of the polynomial, the limiting Newton diagram, 
appears naturally in the context of Renormalized Graeffe Iteration. It allows to recover 
the information about multiple roots and pairs of roots. (See [25]). 

In Section 3, we give a procedure to obtain the limiting Newton diagram of a given 
polynomial. It is effective in the sense that, if we can bound the separation 

101 

max T— 7 , 

ic.i>ioi 101 

then we can effectively identify the multiple roots and pairs of roots. It will converge, 
and eventually provide the list of multiple roots and pairs for any circle-free polyno- 
mial, in a finite (but unknown, not effective) number of iterations. 

The other big weakness of classical Graeffe iteration is the fact that it returns the 
moduli of the roots, but not the actual roots. As a matter of fact, information about the 
argument of the roots is lost, and should be recovered by other means: 

Example 3. Consider polynomials f{x) = a;^ — 2a;-|- 1, g{x) = — 1, h{x) = x'^ + 1. 
After two Graeffe iterations, all the three polynomials are mapped into f{x). 

Many algorithms have been proposed to recover the arguments [28]. In this paper, 
we will differentiate the Graeffe iteration operator, and obtain an iteration defined on 
the appropriate tangent bundle. This new operator will define a mapping between 1-jets 
of polynomials. By the latter we mean expressions of the form f{x) + ef{x), where e 
is a formal parameter. This procedure is discussed in section 4. In the end of the same 
section, we shall discuss the stability properties of this process. 

In section 5, we compare the numerical behavior of Renormjilized Tangent Graeffe 
Iteration to other publicly available algorithms. 



2 Renormalizing Graeffe 

2.1 The Renormalized Graeffe Iteration 

Example 2 shows a typical behavior of classical Graeffe iteration performed by digital 
computers [9] . In order to avoid that sort of overflow, the authors introduced in [22] 
the Renormalized Graeffe Iteration. Although the details and the mathematical foun- 
dations of the algorithm are described in [22J, to keep the present work self-contained, 
we give below a very short description of the main ideas: 

One should consider the computation of g = f as divided in several renormal- 
ization levels. 
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Coefficients of / 
Coefficients of Gf 
Coefficients of G^/ 

Level N Coefficients of f 

At renormaHzation level N, all coefficients gj of g = G^f should be represented 
in coordinates 

rf) =-2-^ log Iff.l, 

and 

"i^^ =9j/\9j\ ■ 

Therefore, we shall obtain convergence of the radial coordinates r^^-* (at least in the 

case of roots of different moduli). The dynamics of the angular coordinates aj^^ is 
typically chaotic. 

In order to pass from level N to level A'^ + 1, a Renonnalized Graejfe Operator was 
defined in [22]. Intermediate computations were performed in coordinates rj^'' and 

a^^^ by means of renormalized arithmetic operations. For instance, the renormalized 
sum (r, a) of (ri , ai ) and (r2 , 0:2 ) can be defined (in renormaUzation level N) by 

r = -2-" log + a^e-^""-' \ 

"~ |aie-2^'-i+a2e-2"'-2| 

Renormalized sum can be computed without overflow by the formula in Algo- 
rithm 1 . This is a simplified, non-optimal version of renormalized sum. Notice that 
one or two of the inputs can be the renormahzation of 0, i.e., cxd. Under the usual con- 
ventions, 00 is greater than any real number. Therefore, if only one of the arguments is 
00, the correct result wiU be returned. 

A few extra mathematical ideas related to the renormalized Graeffe operators, as 
well as some other mathematical results can be found in [22]. 

2.2 The Renonnalized Newton Diagram 

The first goal of this section is to introduce the concept of Renormalized Newton dia- 
gram, which is going to play a key role in the practical implementation of the algorithm 
discussed in this paper. The second is to prove a convergence result based on such idea 
using some earlier results of Ostrowski's. 

We start by reviewing the concept of Newton diagram, which has been used exten- 
sively by Ostrowski, Puiseux and Dumas, among others. 



Level 
Level 1 
Level 2 
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Algorithm 1 RenSum ( ri, ai, r2, a'2> P ) 



{ It is assumed that ri and r2 are real numbers or +00, and that |ai| = \a2\ = 1. 
The number p should be equal to 2^, where N is the renormalization level. This 
routine computes (in renormaUzed coordinates !) the sum of aie~^^^ and a2e~P^^ . } 
if ri = r2 = +00 then 

return +00, 1 
A <— r2 — ri 
if A > then 

t^ai+ a2e-P^ 

returnn- 

else 

t^a2- aieP^ 
retumra- 



Let 



i=0 

be a degree d polynomial. As before, we denote by g = f the iV-th Graeffe iterate 
of/. 

We order the roots of / in nondecreasing order of their moduli, to wit: 

ICi|<---<ICd|. (1) 

If the above inequalities are aU strict, then 



lim 2--log ^ =log(|Ci+i|). 

For each N, consider the piecewise linear function r^^^ : [0,d] — > M U {+00} 
satisfying 

rW(i)^-2-^log|5,|. 

Notice that under the above assumptions, r^^^ is convex for sufficiently large N. (See 
figure 1). Indeed, since |Ci+i| > |Cj|> the inclinations satisfy (for large A/') 

rW(i + l)-rW(i) > rW(i)-rW(i-l). 

It is easy to see that if two consecutive roots, say Q and Ci+i> have approximately 
the same absolute value, then the three corresponding points 

_ 1, r^N) _ 1)^ , ^(A^) (,)^ , + 1, r(^) (, + 1)^ 
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Figure 1: The function r^^^ {i), for iV = 0, 1, 2 



will be approximately aligned. Furthermore, the functions r^^^ converge to a piecewise 
linear convex function. 

However, if the inequalities in (1) are not strict, the functions may fail to 
converge. 

Example 4. Let /(a;) = (a; - l)(x - e'^). Then its N-th Graeffe iterate is 

Therefore, we have r'^^^ (0) = r^^^ (2) = 0, but we also have 

rW(l) = -2-^log|l + e2"'^| = -2-^ log |2 cos 2^-1^1 . 
Depending on the choice of 6, this last value can range anywhere from —2~^ log 2 to 

+0O. 

This is one of the reasons for introducing the convex hull of each r^^\ that will be 
subsequently called the Reno nnalized Newton Diagram. (See figure 2). 

Our approach has the advantage of simphfying some of the arguments by Ostrowski 
in [25] by providing plain convergence of RenormaUzed Newton Diagrams; however, 
we will quote several of the results by Ostrowski in the sequel. 

One of the major goals of Ostrowski in [25] was to obtain effective bounds for the 
moduh of roots. This was possible by introducing of the majorant of a given polyno- 
mial: 

Definition 3. A majorant of a given polynomial is any other polynomial, of same de- 
gree, with normegative coefficients greater than or equal to the given polynomial's 
coefficients. 
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"Convex hull" 
"Original r(l)" 



O 2 4 6 8 10 

Figure 2: The function r^^^ {i) and its Convex Hull 

The first step of Ostrowski's construction is Newton's majorant: 

Definition 4. A polynomial A = X^f^g Aix'^ with nonnegative coefficients is called 
normal if the following conditions hold: 

1. If > and Aj > for « > j, then A; > for alH < i < j. 

2. ¥oil = l,...,d-l, 

Af > Ai-iAi+i . 

A normal majorant T — "^"TiX^ of f is called minimal if for any other majorant T' 
of / we have 

<T;,i = 0,l,...,d (2) 

Notice that Condition 2 above means that the graph of the points of the form 
(/, — log(r;)), for I = 0, . . . , is convex. 
In the language of majorants, 

Proposition 1 (Ostrowski[25]). Any polynomial 

d 
3=0 

possesses a unique minimal normal majorant, 

d 

3=0 
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The polynomial M. / will be called Newton Majorant of the polynomial /. 

The result can be proved by using the convex hull c6 of the function — log \fi\. and 
constructing the polynomial Al/ as the polynomial with positive coefficients {M. f)j = 
e-'!*(j) . We refer the interested reader to Ostrowski's work [25, 26]. 

We remark that if the polynomial has roots of strictly increasing moduli, then the 
coefficients of the Newton's majorant of g = f coincide with \gi\ for N suffi- 
ciently large and i = 0, 1, . . . , d. 

However, the introduction of the Newton Diagram allows us to consider the gen- 
eral situation of possibly many roots of same moduU. As before, we order them in 
nondecreasing order and consider the indices 

io = < ii < i2 < • • • < < = d 

as 0, d and exactly those integers i between 1 and d — 1 such that |Ci| < |Ci+i|- This 
way, we have that 

ICij-i I < ICij-i+l| = • • • = ICij I < ICij+ll • 

The fundamental result, in this case is 

lim 2-^ log = fe+i - I,) log 101 , 

for ij < i < ij^i and < j < I. (c.f. equation (79.8) of [26]). In the language of 
Renormalized Newton Diagrams, that very same equation can be written as: 

log ICil = lim \hl , for u < i < 

N-Kx Ij^i — Ij 

As remarked by Ostrowski, the above formulae are only useful in the determination 
of the moduli of the roots if we know "a priori" the values ii < 12 < ■ ■ ■ < ii- This 
is obviously not the case in most applications. Instead, Ostrowski's results provide 
effective bounds for the convergence of the r^^\ and thus for the values of the | Ci 1' s. 

Theorem 2 (Ostrowski[25],Theorem IX.3). Let 

q{v) ^ 1 - 2-^'^ , 



and 



then 



rp{N) 

rj(N) iJ. ^v-\ 
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As a consequence of the above estimate, Ostrowski gets the following bound 



(2d)-2-" < — — < {2dr 

Corollary 1. If r^^^ (i) denotes the i-th ordinate of the N-th Renormalized Newton 
Diagram, then 



lim frW(i) - rW(i - 1)) = log|Ci 



N~ 

for i = 1, . . . , d Furthermore, the error is bounded from above by 

2"^ log(2(i) . 

This is indeed a strong result, since nothing is assumed on the coefficients or the 
roots of the original polynomial. However, it is possible to get a better error bound by 
assuming a minimal separation on the moduU: 

min M > 1 + e 

iCii>ioi 101 

for some e > 0. 

We note that in the above formula, if e is well defined (i.e. there are at least two 
roots of different modulus) then it is non-zero. 



3 Computation of the Newton Diagram 

3.1 Algorithm and Main Statements 

The main issue in this section is the following: 

We are given a certain polynomial g, obtained after a few Graeffe iterations of a 
polynomial /. The roots of 5 are Zi, ... , Z4, and we order them so that: 

\Zi\ < |^2| < ••• < \Zd\ 

We want to know which of the inequalities are strict. We do not know the actual 
value of the Zj's, we know only the coefficients of g, i.e., the symmetric functions of 
the ZiS. 

We are also ready to assume that 

D • \Zi+l\ 

R = mm , „ , 

\Zi+^\>\Zi\ \Zi\ 

is a large real number. Indeed, if Ci, . . . , Cd are the roots of /, always ordered such as 
ICil < ••• < ICdI.then 



Ki+i|>ICi| ICil 
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is always strictly greater than one. Hence, given any A > 0, by performing N > 
log2 iterations, we can assume that R ^ > A. 

Recall that the RenormaUzed Newton Diagram of g is the convex hull of the func- 
tioni I— > — 2~^log5ij. As A?^ grows, the Renormalized Newton Diagram of 5 converges 
to the convex huU of 

i ^ -log|/o| +^log|Ci| . 

j<i 

However, we want to be able to decide infinite time what are the sharp corners of the 
convex hull of i log \fd\+ J2j<i log I Ci I ) ■ As in the preceding section, we write: 

ri = -2-^log\gi\ 

where g is the N-th iterate of /. Notice that we dropped the superscript N of r\^^ . 

Proposition 2. Let g = f, where f is a degree d polynomial and G denotes the 
Graejfe iteration. Let Z, Q and p be as above. Assume that 

7V>3 + log2— 

logp 

Then, Algorithm 2 below with input N, d, r, p produces the list of the sharp edges of 
the convex hull ofi^ ^j<i 

Remark: Algorithm 2 has running time 0{d). 



3.2 Some Estimates about Symmetric Functions 

In order to prove Proposition 2, we need a few estimates about symmetric functions. 
First of all, let I = {i : \Zi\ < \Zi^i\} U {0. d} be the set of sharp corners of the 
Umiting Renormahzed Newton Diagram. As before, let ak denote the A;-th elementary 
symmetric function, 

'^k{Z)= ^ Zj^...Zj^ 
jl<---<jk 

Then, 

Lemma 1. For iGlwe have 

(Td-i(Z) = Zi+iZi+2 ■ • • Zd{l + c) , 

where |c| < ((^) - l) < 
Proof of Lemma 1: Write 

jl<---<3d-i 
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Algorithm 2 Strict Convex Hull {N,d,r,p) 



{ Create a list A, containing initially the element Aq = 0} 

A,- ^ ; 

{ The error bound below will follow from Lemma 4.} 

E^l (2-^+1 log(2'^ + - 2-^+2 log(l - + i^f^) 

{ Now, we will try to add more points to the list A. At each step, we want to ensure 
that we have always a convex set.} 

for i ^ 1 to d do 

{ We discard all the points in A that are external to the convex huU of A and the 
new point. Let A^ be the last element of A} 



while J > Oand ^^^^3^ > ^73^7 " ^do 

i ^ i - 1 

{ Now, we append the point i} 
Aj ^ i; 
Return (Aq, • • • ,Aj) 



In the sum above, \Zj^ . . . Zj^_.\ < R ^|Zj_|_iZj_|_2 • • • ■^dl for my 
choice of ji, . . . jd-i except i + 1,. ..d. Since there are (^) — 1 other terms, 
we obtain that: 

□ 

Defuiition 5. We wiU say that ii and 12 are successive elements of / if and only if: 

1. ii e / 

2. ii<i<i2=>i^I 

3. 12 e / 

Lemma 2. Let ii and 12 be successive elements of I, and let i\ <l < 12. Then 

\<ja-i{z)\ < (yi^lZ'i) +^') \Zir-'\Zi.+i\---\z<i\ 
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withc'<[{1)-{l-_!l))R-'<2'^R-' 
Proof of Lemma 2: Write 

<^d-i{Z) = ^ Zj^...Zj^_^ 

h<-<3d-i 

= Z^^h--- ^3d-l +2^ Zj^ ... Zj^_, 

where ranges over the j such that ii < ji < • • • < ji2-i < ^2 + 1 and ji^-i+i = 
^2 + 1, . . . , jd-i = d. Of course, ranges over all the other terms. 
We can rewrite as: 

= Zi2+iZi2+2 ■ ■ ■ Zd \ ^ 

\il < jl < • • • < ji2 - i <»2 

Hence, 

The terms in are all smaller than |*^~'|^i2+i| . . . 

Since there are (f) - of them, 

iE"i < ((•) - (t"-/)) 

Adding those two bounds, we obtain indeed: 

|--<^>| ^ ((:-") + ((') - (r-'O) «-') 

□ 

The estimates above can be converted into iogscale' estimates: 

Lemma 3. Let ii, be successive elements of I, and let ii < I < 12. Then the 
following three equations are true: 

1. 

^^^^^^^=log|C.I+2--+Mog(l + c) 

«2 - il 

With \c\ < (ma^ {{t),0) - 1) < 2'^-'- 

2. 

'^'''^ - '-('^ < log IC. I + 2-- log ( h - V) +c']+ 2-- log 11 + c| 



12 — 1 V V ^2 



where \c\ < ((/J - l) R-^ < 2<*i?-i and d < ((^) - {Xll)) ^"^ < 
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3. 

where \c\ < ((/J - l) R'^ < and c' < ((^) - (tl^^)) i?"^ < 

Proof of Lemma 3: By using Lemma 1 with i = 12, we obtain: 

-2-" log Wd-i,{Z)\ = -2-^ (log + ... log |Zd|) - 2-^ log(|l + c"|) (4) 

Using the same lemma with ? = ii, we get: 
-2-^ log \aa-n {Z)\ = -2"^ (log \Zi,+^\ + ... log \Za\) - 2"^ log(|l + c"|) (5) 
Subtracting the two previous expressions and dividing by ^2 — ii we get: 

'^'-'^ '^''^ = 2-^log|Z,J+2-^+ilog(l+c") 
12 - ti 

= logKiJ+2-^+ilog(l + c") 

This shows the first part of the Lemma. 
By using Lemma 2, we can also bound: 

-2-^log|<7d_K^)| > ~2-^{t2-l)\og\Z,,\ 

-2-^(log|Z,,+i| + ...log|Zrf|) (6) 
-2-^log(|(--_7)+c'|) 

where c' is as in Lemma 2. 

We can now estimate equation (4) minus equation (6), altogether divided by 12 — ^: 

< log I + 2-- log ( (^^- + c') + 2-- log |1 + c| 

We can also estimate equation (6) minus equation (5), altogether divided by ^ — ii: 

□ 

3.3 A Decision Criterion 

Lemma 3 can be used do decide if a point in the convex huU of g is converging to a 
sharp corner of the limiting convex hull or not. 

Lemma 4. Assume that 
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a. m> max(i2 — i\) when ii and 12 are successive elements of I. 

b. 2-^+1 log(2™ + - 2-^+2 log(l - < E < 

c. i < j < k 
Then, 

1. Ifi and j are successive elements of I and there is no other element of I between 
j andk, then rUtiriH < ifflzlW _ e 

2. Ifi and k are successive elements of I then > _ ^ 

Proof of Lemma 4: Part 1: Assume that i, j are successive elements of /. Then, part 1 
of Lemma 3 implies: 

''^^\ ~ ""^'^ < log 101 - 2-^+1 log(l - 2'^i?-i) 
3 * 

For the evaluation of ^^^^^^^j^, we have to distinguish two cases: If fc e /, then 

> log lai + 2-^+1 iog(i - 

K — J 

If k /, let m be such that j and m are successive elements of I. Recall that 
j < k < mby hypothesis. Using part 3 of Lemma 3, we get: 

""^^^ ~ "^.^^^ > log ICrr.1 - 2-^ log(2'" + 2'*i?-i) + 2-^ log(l - 2''R-') 
k — 3 

In any case, 

rS^izM < + logic,!- log lai 

J -I k- J 

+2"^ log(2'" + 2'^R-^) - 2-^+2 log(l - 2'^R-^) 

We use the hypothesis E < to deduce that log |Cj| - log |Cfe| + < and: 

^0") - ^ r{k) - r{j) ^ 
j -i k-j 



Part 2: Using Lemma 3, we have: 

^^^^ ~ ^^^^ > log ICfel - 2-^ log(2™ + 2'^R-^) + 2-^ logfl - 2'^R-^) 
3 - « 



r{k) r{j) ^ I ^ ^^^^2" + 2''R-^) - 2"^ log(l - 2''R-^) 



k-j 
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Subtracting, we obtain: 



rjj) - r{i) 
3 - i 



> 



> 



r{k) - r{j) 



- 2 



-AT+l 



k-j 

+2-^+Mog(l - 2''i?-i) 
r{k) - r{j) 



log(2"' + 2''i?-^) 



k-j 



-E 



□ 



3.4 Proof of Proposition 2 

Proof of Proposition 2: Let iV > 3 + loga t^^- It is easy to check that R > 2^<*, 
hence2''i?~^ < 2''^'^. So we can bound: 

2-^+1 log(2™ + - 2-^+2 log(l - 2''i?-i) < 

2(m + 1) logplog2 41og/9 1 logp 
^ 8rflog2 ^ 8rflog2 28-l ^ 

Therefore, we are in the conditions 1 and 2 of Lemma 4, with m = d. Correctness 
of the algorithm 2 can be proved now by induction. 

Induction Hypothesis . At step i, the list A contains Ao,...,As, As+i,Aj where 
Ao, ... As are all successive elements of / and Ag+i, . . . , Aj are not in I. (Possibly, 

we can have s = j). 

The induction hypothesis is true at step 1, with j = 0, and G /. At each step, 
there are two possibilities: 

Case 1: i ^ /. In that case a few of the As_|_i, ■ ■ - Aj may be discarded; but part 1 
of Lemma 4 prevents the algorithm from discarding elements of I. 

Case 2: i e /. In that case, part 2 of Lemma 4 guarantees that all the Ag+i ,.. .Aj 
will be discarded. 

Hence, the induction hypothesis is true at step i + 1. At step d, the last point d is 
added to A. Since de I,A = I. 

A note on the running time: although the usual complexity of a convex hull al- 
gorithm is 0{dlogd) for d points in the plane, the complexity is smaller when those 
points are 'ordered' like ours: {i, r{i)). (Compare with Theorem 4.12 in [30]). Algo- 
rithm 2 has a running time of 0{d) operations (including a fixed number of transcen- 
dental operations). Indeed, each point is added to the Ust A precisely one time. It can 
be discarded only once, so the interior 'while' loop is executed at most d — 1 times in 
one execution of the algorithm. □ 
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4 Tangent Graeffe Iteration 



4.1 Perturbation Methods, Infinitesimals, 1-Jets of Polynomials 

Graeffe iteration provides the absolute values of each root in the case such roots are all 
of different moduli. Recovering the actual value of each root, and recovering pairs of 
conjugate roots or multiple roots require further work. 

Many algorithms have been proposed to recover the actual roots, such as reverse 
Graeffe iteration, splitting algorithms. See [28] and references therein. 

A possibility of theoretical interest would be to consider a perturbation of /; assume 
first that / is a polynomial with roots Ci> • • • > Cd such that |Ci| < IC2I < • • • < |Cd|- 
Then, consider also the iterates of 

fix + e) 

Graeffe iteration of f{x) will provide . . . , K^^l, while Graeffe iteration of 
f{x + e) will provide |Ci — e\, . . . , \Q — e\. Therefore, we will be able to compute: 

ICi|'-|Ci-e|' = -2eRe(Ci)+e2 

thus recovering Q. 

As mentioned before, this is a possibility of theoretical interest only. The per- 
turbation method above would lose half of the working precision in any reasonable 
implementation. Therefore, we will prefer to compute the derivative of |C — with 
respect to e. 

The value e will be treated as an infinitesimal; therefore, instead of storing in mem- 
ory a certain value z + ez, we will store z and z separately. When computing some 
differentiable function G(z + ez), we will obtain a result G{x) + eDG{z)z. So we will 
compute G{z) and DG{z)z, but we will never need to assign an actual value to e. 

A quantity of the form 2; -|- ei is called a 1-jet. It can also be interpreted as an 
element of the tangent bundle of the manifold where z is supposed to live. 

In this paper, we are mainly concerned with 1-jets of polynomials. We will repre- 
sent degree d polynomials as points in M'^+i or C''+^ . Therefore, a 1-Jet of polynomials 
can be represented as a point of M^'^+^ or C^"^"*"^, since we are working with a linear 
space. 

The dot notation (such as in z) will be reserved in this paper to the 'tangent' coor- 
dinate of a 1-jet z + ez. We reserve the notation /' to the derivative ^ / of a univariate 
function / = f{x), and the notation DF to the derivative of a multivariate function 
F. We need the following construction from Calculus on Manifolds [1, 18]: Let G 
be a differentiable function from manifold X into manifold Y. Its tangent map can be 
written, in our 1-Jet notation, as: 

TG: TX ^ TY 

f + ef ^ G{f)+eDGff, 

where as usual TX and TY denote the tangent bundle of X and Y respectively. 

The iteration of the 1-jet of polynomials / -|- ef can be used to recover the actual 
value of the roots of /. For instance: 
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Example 5. Let / be a real circle free polynomial, not vanishing at 0. Consider the 
1-jet f{x + e) = f{x) + e/'(x); its solutions are (j — e, where Q are the roots of /. 
Let g + eg = TG^{f + ef). Lemma 6 below will imply, in the particular case Q is a 
real isolated root, that: 

c,= lim 2--(M_]^ (ii_i2^) 

N^oo \\gj-i\J \gj gj-ij 
In case Q and Cj+i = Cj ^re an isolated pair of conjugate roots, the limit will be: 

ReO= lim ^-N-i(\9^\^ (9^_9J^\ 
\\9j-i\J \9j+i 9j-iJ 

In the following section, we compute the tangent map of the Graeffe operator in 
usual and renormalized coordinates. 

4.2 The Iteration 

Let / + e/ be a 1-Jet of polynomials. Then its Tangent Graeffe Iterate is: 

TG{f{x) + ef{x)) = i-lf (/(v/S) + 6/(VS)) (fi-V^) + ef{-V^)) 

This can be rewritten as: 

TGifix) + efix)) = G{f) + (-l)'^e + 

Precise formulae for computing g + eg = TG{f + ef) are: 

( 9i = + 2E,>l(-l)'+'+V.+,/.-, 

For an efficient root-finding algorithm the equations above need to be renormalized. 
At each step, this is done by replacing products and sums by their renormalized coun- 
terparts. An adjustment is necessary to pass from one renormalization level to another 
(division of the coordinates r by 2). Those adjustments are summarized in Algorithm 3 



4.3 Convergence Results 

It is now time to show convergence of the (Renormalized) Tangent Graeffe Operator. 
Assume one is given a circle free polynomial /. Its roots will be ordered as |Ci| < 
IC2I < • • • Kd|- One can use the (RenormaUzed) Newton diagram to collect together 
the roots with same moduli. Those will represent single roots, multiple roots, or (in the 
case of real polynonoials) pairs of conjugate roots or pairs of multiple conjugate roots. 
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Algorithm 3 TangentGraeffe (N, d, r, a, r, a) 

{ N (Renormalization level) and d (degree) are integers; r and f should be real 
arrays, and a and a should be array of modulus one complex numbers. This rou- 
tine computes, in renormalized coordinates, the Tangent Graeffe Iterate of the 1-jet: 

X^i (ai^^^"^' + eaie"^"'''^ a;'. The coordinates of the result are given in renor- 
malization level N + 1.} 
p ^ 2^+1 
for i ^ to d do 

(r„(-iyaf) 
(si, A) ^ {{n + h)/'2 - log(2)/p, {-iya,a,) 
for j <— 1 to mm{d — do 

{Si,(3i) ^ RenSum(si, f3i, (r^+j + ri_j)/2 + log(2)/p, {-1)''+^ ai+jai^j ,p) 
{si,Pi) ^ RenSum(si,/3i, {n+j + fi-j)/2 + log{2)/p, (-l)'+-'a,+jQ:j_j,p) 
{si,$i) <- RenSum(si,/3i, {ri-j + h+j)/^ + log(2)/p, {-iy-^ai-jai+j,p) 
return (s, (3, s, /3) 



Lemma 5 (Complex case). Let f be a complex circle-free polynomial with roots Ci 7^ 
0, . . . , ordered as in Theorem 1. Let 

*/ • ICi+l| 

p = mm '——^ 
ICi+i|>Kil ICil 

and let g + eg = {TG)^{f + ef). Suppose that j andj + d' are successive elements 
of I ={i:\Ci\< |Ci+i|} U {0; d}. Then, 



Urn --^ f 9j+d' 9j\ _ Cj+d' 



JV^^So d' \gj+d' gjj lO+d'P ' 
Furthermore, the error is bounded by: 

^d+sd ICdl 2",.. ^ 1-1 



1 



Lemma 5 will be proved in Subsection 4.5. 

Also, real polynomials have usually pairs of conjugate roots; they may have pairs 
with multiphcity. In that case, we can show that: 

Lemma 6 (Real case). Let f be a real circle-free polynomial with roots Ci 7^ 0, . . . , 
C,d ordered as in Theorem L Let 

d'f . ICi + l| 



P 



IC.+i|>IC.| IGI 



and let g + eg = {TG)^ f + ef. Suppose that j and j + d' are successive elements of 
I = {i: \Ci\ < |Ci+i|}U{0;rf}. Then, 

lim ^ \ -^"^ ^i+<^' 



N^oo d' \gj+d' gjJ \Cj+d' 
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Furthermore, the error is bounded by: 



The proof of Lemma 6 is also postponed to subsection 4.5. Lemmas 5 and 6 can 
be used to recover the roots of a polynomial from the Tangent Graeffe iterates of its 
1-jet: 



Algorithm 4 RealRecover ( N, d, I, r, a,f,a) 

{ This procedure attempts to recover the roots of the degree d real polynomial 
QfiC"^ '"'a;\ The list of sharp corners of its Newton Diagram is supposed given 
in / = (/o, • • • 1 ^i+size(7))- See Lemma 6 for a justification } 
for A: ^ to Size(/) do 
d' <— h+i - h 

(6, /?) - RenSum (f/,^, - r/,^, , f,, - 2^) 

m ^ exp (2 '"'"+^,"'"''' ) 

X < /3^mexp-2^6 

if Ik+i — h is even and m > \x\'^ then 

y ^ \/m - |a;|2 
else 



for j <- to 7fe+i - /fe - 1 do 
C/fc+j+i =a;+(-lpy 
return ^ 



Algorithm 5 ComplexRecover ( iV, d, /, r, a, f, a ) 

{ This procedure attempts to recover the roots of the degree d complex polynomial 
^ ■ QfiC"^ The list of sharp comers of its Newton Diagram is supposed given 
in /. See Lemma 5 for a justification } 
for /c ^ to Size(/) do 
d! ■*— -ffe+i - Ik 

(6, /?) - RenSum (f,,,, - r,,^„ f,, - -^2^^, ) 

m ^ exp (^2-^^^^^^— 

X < exp— 2^6 

for j <- to 7fe+i - /fe - 1 do 

return ( 
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4.4 The Main Algorithm 

We can now state the algorithm of Theorem 1 . We start with a fixed, arbitrary value for 
= max|f.|>|f^.| . Proposition 2 guarantees that if 

N>3 + log2 ^ 



logp(/) ' 

then after the A'^-th iterate the convex hull of the Newton Diagram of / is computed 
correctly. 

Algorithm 6 Solve (rf, /, isreal) 

{ It is assumed here that / is a degree d, circle-free real or complex polynomial. In 
the general case, one should first find and output the trivial (0 and oo) roots of /, then 
deflate /. After that, one should perform a random real (resp. complex) conformal 
transform on / so it becomes circle-free } 

for i ^ to d do 
if /i ^ then 

^ fi/\fi\ 

else 

ri < log|/i| 

for I <— to d — 1 do 

if// 7^0 then 

else 

aj ■*— 1 
h < — log I //I 
N ^0 
p = 2 
loop 

r, a, f, d ^ TangentGraeffe (A^, d,r,a,f,a) 

I ^ Convex Hull (N, d, r, p) 
if isreal then 

( •*— RealRecover ( A'', d, I, r, a,f,a) 
else 

C ■*— ComplexRecover ( N, d, I, r, a,f,a) 

OutputCi,---Cd 

ifiV>3 + log2 4^then 

{ Proposition 2 implies that at this point, / is indeed correct for all the polyno- 
mials with separation ration > p. Therefore, it is time to decrease p. } 
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After the A^-th iteration, convergence is guaranteed by the following bounds: Ac- 
cording to Lemma 3, at the execution of algorithm Complex Recover (resp. Real Re- 
cover), 



exp ( -. ^-^(54+1 -9iu) 



IWJ(1+'51) 



wherel^il <e2"''+'i°gi+2'^-'". 

Introducing the error bound of Lemma 5, one gets in the complex case: 



N /-. 



with I (^2 1 < 2''+^ W^'^\P ' where a and b are as in the Algorithms. Therefore, 

lOfc+i+i -S| < S\Ci^+j+i\ 

where (5 < {1 + 6if{l + 62) - 1. 

The real case is analogous. According to Lemma 6, 

-N 



9-" Re C- 



so that 



|ReC«,+j+i -x\< (5|CH+i+i| 



and 



where <5' = 5-^0(^2) 

In both cases, \S\ and eventually |^'| are dominated by Since 

the bound in Theorem 1 follows. 

4.5 Proof of Lemmas 5 and 6 

Consider the 1-jet of degree d polynomials / + e f, with solutions Q + cQ. After A'' 
steps of Tangent Graeffe Iteration, we obtain a 1-jet of polynomials 

g + eg={TGf {f + ef) . 

Since the differential of the transformation : ( 1 — > z = C^'^ gives 

DP^ : C + eC ^ {(jf + e2^ {^f ^ , 
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it is clear that the roots + eZj ofg + eg will be 

We can now compute the derivative at e = of + egi = ad-i{Z + eZ). 
Let's denote by Z~ the vector (^i, • • • ,Zj, - ■■ , Z^ e C''"^ 
Take i G / and notice that 



d_ 
de 



Ud-i{Z + eZ) = '^Zjad-i-i{Zq) 

= ^y;ZjCTd-i-i{Z^) 



Thus, 



d_ ad-^{Z + eZ) 
de Od-i{Z) 

e=0 



E 



Zj Zjad-t-i{Zq) 
-' Zj ad-i{Z) 



-gjO-d-»-i(-^) 



Due to Lemma 1 



Zjad-i-i{Z~.) _ sr-Z^ 



3>l 



and 



Z, Z,ad-r-i(,Z.) _ z, Z^ 

^z, ad-i{z) - Z^- - k^ + m)^ 



j<i 



j<i 



Z-i Z,. 



j 



where 



\vj\<^Qr-'. 



Although this estimate is by no means sharp, it suffices for our purposes. 
Using that i G I and hence \zj/zi+i\ < for i < j, we get 



d_ cTd-i{z + ez) 
de (Td-i{z) 



< max 



R 



< max 
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Therefore, if we take the logarithmic derivative of the expression 
(5 + e5) = (TG)^(/ + e/) 
and evaluate at a successive couple of elements i\,i2 € / we get 



( ih. 




\9ti 


9iJ 



il<j<i2 



< max 

R I 



Now let's assume we are under the hypothesis of Lemma 5. Then, since ii and i2 
are consecutive, it follows that 



ii<j<i2 ii<j<i2 



and so 



-N 



9i2 9ii 



i2 - ii \9i2 9n 



I Ci2 I Ci2 



\Ci2 



-1 2'^+3ci Id 
-rt(^2-^l) ICsl 



(7) 



On the other hand, if we are under the hypothesis of Lermna 6, we get 

h<j<i2 l''*^! 

Thus, 



-N 



i2 - ii \g 



9ii 



■ReC, 



10: 



< 'm'- ^ T7 

i^(^2 - ii) '■.s ic 



^ • (8) 



In either case, we have that each right hand side of equations (7) and (8) is bounded 



by 



d 



R{i2-ii) ICil ■ 
Now, using the fact that the separation radius at the A^-th step 

R = P^ , 

where p is the separation radius of the original roots of / before applying Graeffe, we 
get that 



R{i2 - ii) ICil 



«2 - k ICil 



, as iV 



□ 
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4.6 'Deterioration of Condition' and Stability Properties 

It is important to understand that we will never have to solve g{x) = {G^ f){x). 
Therefore, the actual condition number of g does not matter at all. In order to de- 
termine the roots of /, we will be using the extra information provided by g, where 
{g,g)=TG^{f,f). 

A vaUd source of concern is the propagation of rounding-off error. In the Tangent 
Graeffe algorithm, that error would typically double at each step (it actually doubles at 
each step 'in the limit'). 

However, Lemmas 5 and 6 guarantee that the truncation error decreases as p . 
Hence, in order to obtain a truncation error smaller than a certain 5 > 0, we need 
N X c + log2 log2 S~^, where c = — log2 log2 p is a constant depending only on the 
original polynomial /. 

In order to reduce the accumulated rounding-off error to the same order, one would 
need that 

Ce^2^ < 6 

where is the 'machine epsUon' and C is a constant depending on /. Thus, we just 
need 

e <i2-^« - 

em<^^ ~ C2'=log2(5-i 

What is a reasonable value for S ? The strength of Renormalized Tangent Graeffe 
Iteration is its capacity to solve the 'global' problem: given a polynomial, approximate 
all its roots. Once a suitable approximation of each root is found, local iterative algo- 
rithms (such as Newton Iteration) cheaply provide better refinements of the roots. Such 
a two-step procedure entails a reasonable range of values for 6. Namely, 6 should be 
smaller (but not much smaller) than the radius of quadratic convergence of Newton's 
iteration. 

This radius is of the order of the reciprocal condition number of the original poly- 
nomial. (See [3] Theorem 1 and Remark 1 p. 263 for a precise statement). 



5 Numerical Results and Final Remarks 

5.1 Numerical Results 

A polynomial solver based on the algorithms above was implemented and tested un- 
der IEEE 754 double and double-extended arithmetic. The results below are intended 
to make a case in favor of the stability and the practical feasibility of Renormalized 
Tangent Graeffe Iteration. 

The first set of tests was designed to measure the performance of our algorithm for 
large degree polynomials. The test polynomials are pseudo-random real (Table 1 and 
Figure 3) and complex (Table 2 and Figure 4) polynomials, under the U (2)-invariant 
probability measure [17, 34, 22]. Under this probability measure, random polynomials 
are well-conditioned on the average. 
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"Renormalized Tangent Graeffe" 
"Jenkins and Traub" 



Figure 3: Real pseudo-random polynomials 

The results were certified using alpha- theory [38, 33, 20]. The run- 
ning time (certification excluded) was compared to the code of Jenkins and 
Traub [14, 13] for the values where this code succeeds. 

Running time is measured in user- time seconds of a Pentium- 133 computer miming 
Linux and the gcc compiler. 

In the second set of experiments, we tried to check the behavior of Renormalized 
Tangent Graeffe Iteration in the presence of very badly conditioned polynomials. The 
test polynomials are Wilkinson' s 'perfidious' polynomials [41] 



Pd{x) = {x - l){x - 2) ■ 
and Chebyshev polynomials (Table 3). 



{x - d) 



0<m<d 



2d d 



m) 



The error of the solutions of the perfidious polynomials is measured as max | C, — 
round Cl^ C ^ 'solution' found by the program. Similarly, the error in Chebyshev poly- 
nomials is measured as max \m— round m\ where m = '^'^''^^"^^ ? ^^^j ^ ^ 'solution' 
found by the program. Again, those results are compared to the ones provided by the 
software by Jenkins and Traub. 



5.2 Further Practical Remarks 

• Graeffe process (and hence our algorithm) is known to be parallelizable [12,31]. 

• The algorithm presented here needs only 0{d) memory storage; therefore, all 
intermediate computations for a reasonable degree d fit into the 'cache' memory 
of modem computers. 
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Degree 


Algorithm 





1 




'} 


4 


Seed 

5 


6 


7 


8 


9 
























0.07 




J-1 


0.03 


0.01 


0.04 


0.02 


0.0.^ 


0.0.^ 


0.0/ 


0.04 


0.04 


0.03 


100 


Graeffe 


0.20 


0.20 


0.20 


0.21 


0.21 


0.21 


0.19 


0.20 


0.19 


0.20 




J-T 


0.12 


0.11 


0.12 


0.09 


0.12 


0.12 


0.14 


0.09 


0.10 


0.10 


150 


Graeffe 


0.36 


0.36 


0.36 


0.36 


0.37 


0.35 


0.36 


0.37 


0.35 


0.36 




J-T 


0.26 


0.26 


0.28 


0.24 


0.23 


0.24 


0.27 


0.20 


0.26 


0.24 


200 


Graeffe 


0.56 


0.54 


0.56 


0.56 


0.55 


0.55 


0.56 


0.55 


0.56 


0.55 




J-T 


0.53 


0.42 


0.51 


0.35 


0.43 


0.38 


0.54 


0.40 


0.47 


0.36 


250 


Graeffe 


0.77 


0.78 


0.77 


0.78 


0.78 


0.78 


0.78 


0.78 


0.77 


0.76 


300 




1.02 


1.02 


I.OI 


f.02 


1.02 


1.00 


1.02 


1.01 


1.01 


1.03 


350 




1.29 


1.29 


1.29 


f.29 


1.30 


1.29 


1.29 


1.29 


1.28 


1.29 


400 


Graeffe 


1.59 


1.59 


1.57 


1.58 


1.58 


1.60 


1.57 


1.59 


1.59 


1.58 


450 


Graeffe 


1.90 


1.90 


1.89 


1.91 


1.89 


1.89 


1.90 


1.90 


1.91 


1.91 


500 




2.32 


2.27 


2.31 


2.30 


2.31 


2.31 


2.30 


2.29 


2.32 


2.34 


550 


Graeffe 


2.57 


2.58 


2.58 


2.58 


2.59 


2.59 


2.58 


2.59 


2.57 


2.59 


600 


Graeffe 


2.96 


2.98 


2.97 


2.96 


2.95 


2.97 


2.96 


2.97 


2.98 


2.96 


650 


Graeffe 


3.38 


3.39 


3.37 


3.36 


3.37 


3.34 


3.38 


3.37 


3.37 


3.37 


700 


Graeffe 


3.78 


3.78 


3.78 


3.80 


3.79 


3.77 


3.78 


3.80 


3.78 


3.77 


750 


Graeffe 


4.21 


4.22 


4.21 


4.21 


4.19 


4.20 


4.19 


4.22 


4.20 


4.21 


800 


Graeffe 


4.66 


4.66 


4.65 


4.65 


4.64 


4.64 


4.63 


4.65 


4.65 


4.65 


850 


Graeffe 


5.18 


5.20 


5.18 


5.17 


5.23 


5.18 


5.19 


5.20 


5.19 


5.17 


900 


Graeffe 


5.61 


5.61 


5.58 


5.59 


5.60 


5.57 


5.56 


5.60 


5.60 


5.60 


950 


Graeffe 


6.16 


6.16 


6.15 


6.16 


6.13 


6.14 


6.16 


6.18 


6.14 


6.18 


1000 


Graeffe 


6.60 


6.58 


6.58 


6.58 


6.59 


6.61 


6.59 


6.60 


6.58 


6.59 



Table 1: Real pseudo-random polynomials 



"Renormalized Tangent Graeffe" 
"Jenkins and Traub" 



Figure 4: Complex pseudo-random polynomials 
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Running time (s) 



Degree 


Algorilhm 










Seed 



















1 


2 


3 


4 


5 


6 


7 


8 


9 


50 




0.13 


0.12 


0.12 


0.13 


0.12 


0.10 


0.11 


0.10 


0.11 


0.11 




J-T 


0.07 


0.08 


0.09 


0.09 


0.07 


0.08 


0.07 


0.08 


0.07 


0.09 


100 




0.34 


0.32 


0.32 


0.32 


0.33 


0.33 


0.30 


0.33 


0.31 


0.32 




J-T 


0.25 


0.28 


0.26 


0.25 


0.28 


0.26 


0.26 


0.27 


0.26 


0.25 


150 




0.60 


0.60 


0.61 


0.60 


0.61 


0.60 


0.60 


0.61 


0.60 


0.60 




J-T 


0.55 


0.60 


0.58 


0.66 


0.61 


0.60 


0.60 


0.59 


0.59 


0.59 


200 


Graefft; 


0.95 


0.95 


0.97 


0.94 


0.96 


0.94 


0.92 


0.94 


0.94 


0.93 




J-T 


1.13 


1.17 


1.07 


1.04 


1.09 


1.06 


1.05 


1.06 


1.03 


1.05 


250 




1.33 


1.34 


1.30 


1.32 


1.35 


1.32 


1.32 


1.32 


1.35 


1.33 




J-T 


1.63 


1.68 


1.64 


1.68 


1.66 


1.70 


1.66 


1.63 


1.63 


1.66 


300 




1.77 


1.76 


1.77 


1.78 


1.77 


1.76 


1.77 


1.75 


1.79 


1.76 




J-T 


2.34 


2.49 


2.43 


2.44 


2.50 


2.46 


2.48 


2.45 


2.36 


2.51 


350 


Graefft; 


2.26 


2.30 


2.26 


2.25 


2.28 


2.34 


2.08 


2.36 


2.27 


2.26 




J-T 


3.30 


3.41 


3.39 


3.34 


3.36 


3.72 


3.51 


3.45 


3.35 


3.52 


400 




2.74 


2.75 


2.77 


2.76 


2.77 


2.76 


2.75 


2.77 


2.74 


2.74 


450 




3.27 


3.28 


3.28 


3.30 


3.38 


3.30 


3.30 


3.28 


3.34 


3.30 


500 




3.92 


3.87 


3.89 


3.89 


3.88 


3.89 


3.91 


3.87 


3.90 


3.91 


550 




4.49 


4.48 


4.51 


4.52 


4.50 


4.52 


4.51 


4.49 


4.50 


4.48 


600 




5.16 


5.14 


5.17 


5.18 


5.15 


5.14 


5.15 


5.15 


5.18 


5.13 


650 




5.89 


5,88 


5.87 


5.83 


5.83 


5.86 


5.89 


5.84 


5.89 


5,86 


700 


Graefft; 


6.62 


6.59 


6.59 


6.58 


6.59 


6,61 


6.59 


6.59 


6.58 


6,62 


750 




7.37 


7,42 


7.40 


7.32 


7.48 


7.27 


7.30 


7.30 


7.43 


7.31 


800 


Graefft; 


8.10 


8.06 


8.07 


8.03 


8.10 


8.07 


8.10 


8.06 


8.11 


8.09 


850 




8.96 


8.94 


8.92 


8.95 


8.97 


8.94 


8.92 


8.92 


8.95 


8.96 


900 


Graefft; 


9.72 


9.70 


9.70 


9.70 


9.72 


9.72 


9.71 


9.70 


9.67 


9.71 


950 


Graeife 


10.63 


10.67 


10.65 


10.64 


10.67 


10.65 


10.67 


10.61 


10.65 


10.60 


1000 


Graeffe 


11.49 


11.47 


11.54 


11.46 


11.50 


11.49 


11.54 


11.50 


11.51 


11.50 



Table 2: Complex pseudo-random polynomials 



Perfidious polynomials 



Degree 


Algorilhm 




10 


Grat;ffe 
J-T 


5.123013 X 10"-^^ 
4.859935 X 10"^^ 


15 


Grat;ffe 
J-T 


3.968295 X 10~*'^ 
5.508868 X 10^"^ 


20 


Grat;ffe 
J-T 


1.780775 X 10~*'-^ 
1.275754 X 10~04 



Chebyshev polynomials 



Dt;groo 


Algorilhm 


Error 


10 


Graefft; 
J-T 


8.790711 X 10"-^*^ 
8.790711 X 10~^*^ 


15 


Graeffe 
J-T 


2.169163 X 10"-^'"^ 
2,169163 X in~^"' 




J-l 


1 9(,:iM,. y lu 
i,27.s977 X 10 


25 


Graeffe 
J-T 


1.266375 X lO"-'--'- 
1.663025 X 10~11 


30 


Graeffe 
J-T 


5.511325 X 10~^^ 
3.301608 X 10"^*^ 


35 


Graeffe 
J-T 


5.708941 X 10~^y 
3.840000 X 10"*^^ 



Table 3: Perfidious and Chebyshev polynomials 
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• Use of higher precision may be required to handle very badly conditioned poly- 
nomials; as a matter of fact, polynomials of that sort are only meaningful if their 
coefficients are known with sufficiently high accuracy. For example, when they 
are obtained by symbolic manipulation. 
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